Time-dependent Gutzwiller theory of magnetic excitations in the Hubbard model 
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We use a spin-rotational invariant Gutzwiller energy functional to compute random-phase- 
approximation-like (RPA) fluctuations on top of the Gutzwiller approximation (GA). The method 
can be viewed as an extension of the previously developed GA+RPA approach for the charge sector 
[G. Seibold and J. Lorenzana, Phys. Rev. Lett. 86, 2605 (2001)] with respect to the inclusion of 
the magnetic excitations. Unlike the charge case, no assumptions about the time evolution of the 
double occupancy are needed in this case. Interestingly, in a spin-rotational invariant system, we 
find the correct degeneracy between triplet excitations, showing the consistency of both computa- 
tions. Since no restrictions are imposed on the symmetry of the underlying saddle-point solution, 
our approach is suitable for the evaluation of the magnetic susceptibility and dynamical structure 
factor in strongly correlated inhomogeneous systems. We present a detailed study of the quality of 
our approach by comparing with exact diagonalization results and show its much higher accuracy 
compared to the conventional Hartree-Fock+RPA theory. In infinite dimensions, where the GA be- 
comes exact for the Gutzwiller variational energy, we evaluate ferromagnetic and antiferromagnetic 
instabilities from the transverse magnetic susceptibility. The resulting phase diagram is in complete 
agreement with previous variational computations. 



PACS numbers: 71.10.-w, 71.27. +a, 71.45.Gm 

I. INTRODUCTION 

It is now about 40 years ago that Gutzwiller proposed a 
variational wave function for correlated electronic models 
with a purely local interaction, i.e., for the Hubbard- like 
modelsji*£ The basic idea is to partially project out con- 
figurations with doubly-occupied sites from the Fermi sea 
in order to optimize the contributions from kinetic and 
potential energy. As a consequence, in contrast to the 
conventional Hartree-Fock (HF) theory, the Gutzwiller 
wave function captures correlation effects like the band 
narrowing already on the variational level. However, the 
exact evaluation of the ground state energy within the 
Gutzwiller wave function is fairly difficult and up to now 
has only been achieved in one and infinite dimensions^ In 
the latter case the solution is equivalent to the so-called 
Gutzwiller approximation (GA) which has been applied 
to describe a variety of finite dimensional systems rang- 
ing from the properties of normal 3 He (see Ref. 0) to the 
stripe phase of high-T c cupratesA^ 

The GA in its original formulation was restricted to ho- 
mogeneous paramagnetic systems and only later on gen- 
eralized to arbitrary Slater determinants by Gebhard 7 
and, more recently, by Attaccalite and Fabrizioi& The 
same energy functional was obtained from the Kotliar- 
Ruckenstein (KR) slave-boson formulation of the Hub- 
bard model when the bosons are replaced by their mean- 
values^ Moreover, the KR slave-boson approach provides 
a controlled scheme for including fluctuations beyond the 
mean-field solution. Formally this has been achieved 
by several authors within the functional integral formal- 



j sm> io i ii i i2 J i3 1 i4 u oweverj the expansion of the KR hop- 
ping factor z SB turned out to be a highly nontrivial task, 
both with respect to the proper normal ordering of the 
bosons and with respect to the correct continuum limit of 
the functional integral,— These difficulties have severely 
hampered the computation of charge fluctuations within 
the slave boson approach. To our knowledge this tech- 
nique has therefore only been applied to toy models 1 ® 
and to compute the optical conductivity in the param- 
agnetic regime j 13 : 14 The latter, however, did not lead to 
controlled sum rules due to the above mentioned difficul- 
ties iA£ 

In Refs. H^IT^ we have developed an alternative scheme 
for the computation of random-phase-approximation-like 
(RPA) fluctuations beyond the GA. Our approach, la- 
beled GA+RPA, is based on well developed techniques 
in nuclear physics 1 ^ and RPA fluctuations are obtained 
in the small oscillation limit of a time-dependent GA. By 
comparing with exact diagonalization results, we have 
shown that the computation of static and dynamical 
correlation functions performs much better within the 
GA+RPA than within conventional HF+RPA theoryiS 
Since no restrictions are imposed on the symmetry of 
the saddle-point solution, the GA+RPA method is also 
suitable for the investigation of strongly correlated elec- 
tronically inhomogeneous systems. Based on this formal- 
ism, two of us have recently explained the evolution of 
the optical conductivity with doping in high-T c cuprate 
compounds i2£ 

Our previous investigations were restricted to the eval- 
uation of RPA fluctuations in the charge sector where 
the total spin is conserved by the particle-hole excita- 
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tionsi«iii However, in general, one has to distinguish 
between longitudinal (i.e., with AS Z — 0) and transverse 
spin excitations (i.e., with A5 2 = ±1), the latter involv- 
ing particle-hole pairs with opposite spins. Therefore, 
longitudinal excitations are optically allowed by dipole 
selection rules whereas transverse excitations can be ex- 
cited by spin-carrying particles like neutrons. For spin- 
rotational invariant systems, the triplet transverse exci- 
tations with AS Z = ±1 are degenerate with the triplet 
longitudinal excitation with AS Z = and, therefore, it is 
enough to solve the problem in the longitudinal channel. 
As discussed below, the solution in the transverse channel 
is useful for a consistency check. If spin-rotational sym- 
metry is broken, (e.g., for ferromagnetic or spin-density- 
wave states) the triplet excitations will split and one has 
to solve both channels to obtain the whole spectrum. 

The present paper is therefore devoted to the computa- 
tion of transverse magnetic excitations on top of the GA. 
Various approaches have been already adopted in order 
to accomplish this task. In Ref.0, Biinemann has evalu- 
ated the spin-wave excitations in itinerant ferromagnets 
by determining variationally the energy of the excited 
state S+I^g)) where \^g) denotes the Gutzwiller wave 
function and is the spin-flip operator with momentum 
q. Furthermore, spin excitations around para magnetic 
saddle points have been investigated in Refs. Illll2l22l 
within the functional integral technique based on the 
spin-rotational invariant slave-boson schemed 

Our investigations below are related to these previous 
investigations but differ in two important aspects. First 
we will eliminate the bosonic degrees of freedom (except 
for the double occupancy D) from the energy functional, 
which thus only depends on the density matrix and the 
parameters D. Formally, this procedure defines an ef- 
fective Gutzwiller Hamiltonian, which can be expanded 
with respect to both charge and spin fluctuations. As 
usual, both types of excitations are decoupled in case of 
saddle points with collinear spin structure. Second, the 
density matrix can be constructed from arbitrary Slater 
determinants, and, therefore, the method is suitable for 
the investigation of magnetic excitations in inhomoge- 
neous systems. In this respect, the size limitations in 
numerical solutions are exactly the same than for the in- 
homogeneous HF+RPA approach^ 

The paper is organized as follows. In Sec. [H] we de- 
rive the GA energy functional from the spin-rotational 
invariant slave-boson Hamiltonian and show how RPA 
fluctuations in the charge and spin channel can be ob- 
tained within the time-dependent Gutzwiller approach. 
In particular, we focus on the magnetic excitation spec- 
trum obtained in this way from the Hubbard model. Re- 
sults for specific systems are presented in Sec. IIIII As 
a first example, we consider in Sec. IIII Al the two-site 
Hubbard model, where the analytical solution is avail- 
able for comparison. Since at small U the mean field 
ground state is spin-rotationally invariant, the expected 
degeneracy between longitudinal and transverse spin ex- 
citation allows us to check the consistency among charge 



and magnetic channel computations. Then, in Sec. IIII 51 
the method is applied to a homogeneous and paramag- 
netic GA solution, where it turns out that the evaluation 
of transverse magnetic susceptibilities is greatly simpli- 
fied as compared to previous approaches. In particular, 
we evaluate the ferromagnetic and antiferromagnetic in- 
stability lines for an infinite-dimensional hypercubic sys- 
tem, and demonstrate the exact agreement with varia- 
tional results. Section ITlI CI is devoted to a comparison 
of the GA+RPA magnetic excitation spectra with exact 
diagonalization and HF+RPA results respectively. Con- 
cluding remarks appear in Sec. IIVI 

II. FORMALISM 
A. Spin-rotational invariant GA 

The starting point is the one-band Hubbard model: 

H = Y^ %4,a c j> + U H n *.T n a> (!) 

where Ci^ a (cJ CT ) destroys (creates) an electron with spin 

a at site i, and 7i; iCr = cl a Ci y(T . U is the on-site Hub- 
bard repulsion and tij denotes the hopping parameter 
between sites i and j. Our investigations are based on 
the spin-rotational invariant form of the slave-boson ap- 
proach introduced by KRm& Within this formalism one 
introduces auxiliary bosons (e\) and di (d\) which rep- 
resent the annihilation (creation) of empty and doubly- 
occupied sites, respectively. In addition, the single oc- 
cupied states are represented by two particles, a spin- 
1/2 fermion and a boson p which can have either spin 
S = or S = 1 in such a way that the combination 
has spin-1/2. The four p states (a singlet and a triplet) 
are combinations of the elements pi^a 1 of a 2 x 2 matrix 
Pi. In the saddle-point approximation all boson opera- 
tors are treated as numbers and the matrix can be 
parametrized as: 

( Pa ^exp(-#,)\ 
V 7f^ ex P(+*^) ) ' 

with pi, p itT , and 4>i real. 

Besides the completeness condition 

e t 2 + tr(p* Pi )+A = l, (3) 
the boson fields are constrained by the following relations 

tr(vp? Pi ) + 2^,0 A = ^(r^p'f', (4) 

where, in general, = {c\ a Cj, a i) denotes the density 
matrix, are the Pauli matrices (including r = 1), and 
A = dl 
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After rewriting the Hamiltonian Q in terms of fermion 
and boson operators^ we can construct a spin-rotational 
invariant Gutzwiller functional by eliminating the boson 
fields except for A via the constraints, Eqs. J3J and 10}. 
As a result, one obtains: 



e ga _ \- , 
where the matrix z, reads as 



CT1,(T2 



C/^A, (5) 



Within this notation we can formally write the GA en- 
ergy as 



(13) 



where k = p,h labels particle and hole states and tk are 
the corresponding one-particle energies. 



Zj + cos 2 j + Zi sin 2 -^r [zi + — z t ] cos $ 



s 



'- [zi + — Zi ] cos $ Zi + sin 2 ^ + Zi cos 2 $ 



with 
tan 2 $ 



±\2 



(Af) 



^(l-A-lAf) 2 ) (^-A-(A?) 2 ) 
- A ± 5fVl + tan 2 $, 



,(8) 
(9) 



and for clarity spin expectation values are denoted by 
Sf = Pli 1 , Sr = p±;\ S? = (p]f - P j; l )/2, and 
pa = p\i + pli 1 . Note that in the limit Sf = 0, 
where the matrix z^ is diagonal, one recovers the stan- 
dard Gutzwiller energy functional as derived by Geb- 
hardi or KRiS Furthermore, it has been shown that the 
spin-rotational invariant slave-boson scheme can be de- 
rived from the KR (or alternatively Gebhard's) energy 
functional when the spin rotation is applied to the un- 
derlying Slater determinant^ Therefore, Eq. (|3J) can be 
viewed as the more general GA-like energy functional for 
a Hubbard Hamiltonian. 

In order to obtain the stationary solution of Eq. © 
one has to minimize E GA with respect to the double oc- 
cupancy parameters D and the density matrix p. The 
latter variation has to be constrained to the subspace of 
Slater determinants by imposing the condition 



P = P: 



(10) 



which is equivalent to the diagonalization of the elec- 
tronic problem supplemented by the variation with re- 
spect to D only. A detailed description of the corre- 
sponding formalism can be found in Ref. I2M 

Regarding the stationary solutions, we will restrict to 
Slater determinants which are diagonal in spin space, 

i.e., plf ^ = pij 'S^o-t . Thus we do not consider spin 
canted solutions 2 - which would mix charge and spin ex- 
citations. Therefore, the diagonalized density matrices 
have eigenvalue 1 below the Fermi level (= hole states: 
h) and zero above (= particle states: p) and consequently 
are also diagonal in spin space: 



(0) 
Pha,hc 

«(°) 



(11) 

(12) 



B. Calculation of RPA fluctuations around general 
GA saddle points and magnetic excitations 

The energy functional Eq. JSJ) is a convenient starting 
point for the calculation of charge and spin e xcitati ons 
on top of general GA wave functions. In Refs. Il6ll7l we 
have already given a detailed derivation of the GA+RPA 
formalism in the charge sector, which, in the following, 
we extend to include the spin fluctuations. 

We thus study the response of the system to an exter- 
nal time-dependent perturbation 

F(t) - E [k^'iM^+H.c], (14) 

i,j,<r,a' 



fij,<r<r f (t) — fij.crcr' (0)^ 



(15) 



which induces small amplitude oscillations of D and p 
around the GA saddle point: 



D = D {a) +5D(t), 
P = pW+6p(t). 



(16) 
(17) 



Correspondingly, we have to expand the energy func- 
tional Eq. (JSJ) around the stationary solution up to second 
order in the density- and double-occupancy deviations. 
Due to the fact that we restrict to collinear saddle-point 
solutions, the charge and spin sectors in the expansion 
are decoupled and one obtains 



E[p, D}= E + ti{h°Sp} + SE 



charge 



SE spv 



(18) 



where the subscript indicates quantities evaluated 
in the stationary state and we have introduced the 
Gutzwiller Hamiltoniani 27 ! 28 



Kf [p,D} = 



dE GA 

; ( 

dpa 



(19) 



SE charge contains the expansion with respect to the 
double-occupancy parameters and the part of the density 
matrix, which is diagonal in the spin indices. This part 
of the RPA pr oblem has already been studied in detail 
111 Refs. Il6ll7l where it was shown that the SD fluctu- 
ations can be eliminated by assuming that they adjust 
instantaneously to the evolution of the density matrix 
(antiadiabaticity condition). 
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The spin part of the expansion reads: 

id><? 



(20) 



hi," 



with the following abbreviations for the quadratic parts 
of the z-factor expansion 



9pl 



'-bp. 



a. (J rit ' 



ft 1 ? 

U ^1,(7,(7 



a.— a c —a.tj 



> c a, — o r 

=^ s Pu °P 



(21) 



(22) 



The explicit results for the derivatives are given in ap- 
pendix^ It is interesting to observe that, in contrast 
to the charge excitations, the evaluation of the magnetic 
excitations can be performed without any adjustment of 
8D to 5/9, i.e., without any assumption on the time evolu- 
tion of D. Only in the case of non-collinear saddle points 
one would have a coupling between spin and charge fluc- 
tuations and, therefore, the necessity to invoke the an- 
tiadiabaticity condition to eliminate the SD deviations. 

The density fluctuations Sp in the expansion Eq. i|18|) 
are restricted to the subspace of Slater determinants, i.e., 
they have to obey the constraint Eq. (|10|) . One can 
therefore divide Sp into the particle (p) and hole (h) sec- 
tors using the property of the density matrices Eqs. i|TT|) 
and fT5]b: 



{sp h £} 



„(o) 



(i-pB)Spp { Z, 

(1 



P { °1)6 ■ • • 



(23) 
(24) 
(25) 
(26) 



where by {Sp^.,} we mean a matrix whose non-zero 
generic elements are of the form Sp^j,. Moreover, one 
can show [see Eqs. (34)-(36) in Ref. Il7j that the pp and 
hh density projections yield a quadratic contribution in 
the ph and hp matrix elements in the small amplitude 
approximation: 



6f& « -J25p h Z„8f?Xn 

pa" 

6pZ, « ^6fC,,5pX>- 



(27) 
(28) 



ha'' 



Hence, although the Gutzwiller Hamiltonian Eq. I|19|) 
is diagonal in spin space, it turns out that the term 
ti(h°dp) — Y^n ^vPw m Eq. (fT%)l (which is first order 
in the pp and hh density projections) yields a quadratic 



contribution in the ph and hp matrix elements: 



hh 

J aai 



tr(h°Sp) = "£e p SpZ+J2e h Sp h a 

pa ha 

= J2(e p -e h )SfC,6p%. (29) 

phaa' 

The fluctuations which are diagonal in the spin indices 
{8p ph a and Sp^ p ) contribute to the expansion in the charge 
channel^ whereas the non-diagonal elements describe 
the zero-order (non-interacting) spin-flip excitations of 
the saddle-point Slater determinant. 

Thus, up to second order in the particle-hole (spin) 
density fluctuations, one obtains for the energy expan- 



5E sp ln 



-(Sp hp ,Sp ph ) 



A 

B* 



B 

A* 



SpP' h ' 
Sp h ' p ' 



(30) 



where the explicit expressions for the RPA matrices A 
and B are given in the appendix [BJ Note that the 
shorthand notation in Eq. (|30l) and below implies that 
p- and h-states have opposite spin, i.e., 5p ph represent 



the joint set of elements of types Sp^ and Sp p 'l- Fol 



O" TJ. '4-1 
we can now evaluate the response func- 
tion corresponding to the perturbation Eq. i|14fl . In case 
of non-diagonal perturbations (as the coupling to a cur- 
rent), one has to define an associated Gutzwiller operator 
which contains the GA hopping matrices Zj. However, in 
the spin channel the most relevant perturbations couple 
an external field locally to some spin operator. The field 
fij,aa' — fii.a-a-'Sij is therefore diagonal in the site repre- 
sentation and remains unchanged within the GA. Upon 
transforming the perturbation to the particle-hole rep- 
resentation one can derive the following linear response 
equation: 



'IV 



A 
B* 



B 
1 



— h\o 



o 



Sp ph 
6p hp 



fph 
fhp 



The inversion of Eq. (|31|l yields a linear relation between 
the external field and the change in the density: 



6p = R(u)f, 



(32) 



and defines the linear response function R(oj) which in 
the Lehmann representation reads as: 



R[V)ph,p'h' 



E 

n>0 



^ph^p'h' 

uj — Q n + ie ijj + fl n + ie 



1 p'h' 1 ph 



(33) 



where we have introduced the eigenvectors of the RPA 
matrix: 



(0|o^fltp|n) 
(0\a p a h \n) 



^phi 
I hp- 



(34) 
(35) 



and |n) denote the unprojected (i.e., without Gutzwiller 
correlations) excited states of the RPA problem. 
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III. RESULTS 

The RPA formalism derived in the previous section 
constitutes a convenient starting point for the calcula- 
tion of spin excitations on top of the GA. One of the 
advantages of the present approach is that it is suitable 
for general Slater determinants, i.e., without any restric- 
tion on translational and (longitudinal) spin symmetries. 
The system sizes which can be treated are the same than 
for the traditional HF+RPA approximation. 19 However, 
also for homogeneous and paramagnetic saddle points the 
GA based RPA approach provides a convenient method 
for the evaluation of spin fluctuations. Our method is 
solely based on the expansion of the density matrix in 
terms of particle-hole fluctuations and does not involve 
other degrees of freedom as in the related functional in- 
tegral slave-boson schemeiii*i2i2£ First, this advantage is 
demonstrated for a two-site Hubbard model which is also 
a convenient toy model for the RPA formalism derived in 
the previous section. However, the GA+RPA approach 
can also be applied within the more conventional Green's 
function technique which is used in Sec. 1111 51 to evaluate 
spin susceptibilities for a homogeneous and paramagnetic 
hypercubic lattice in infinite dimensions. In this case the 
GA becomes exact for the energy functional within the 
Gutzwiller wave function, and we recover the magnetic 
instability lines determined previously by Fazekas and 
collaborators^ The remainder of this section is then de- 
voted to a detailed analysis of the quality of our approach 
by comparing with HF+RPA and exact results for small 
clusters, where the exact solution is known by exact di- 
agonalization techniques. 



A. Two-site Hubbard model 

As a first example, we consider the two-site Hubbard 
model at half filling which can be solved exactly and 
can be studied analytically with both the GA+RPA and 
HF+RPA approximations. On general grounds a mean- 
field (or time-dependent mean-field) approach is expected 
to improve as the dimensionality of the space increases, 
and, therefore, this zero-dimensional problem is the worst 
case and may give an estimation of the maximum error 
which can be expected for these mean- field approaches. 

The exact ground-state energy is given by: 



U - \J~U 2 + 16i 2 



(36) 



and the corresponding eigenfunction reads as 
l*o) - a-j= [| Ti-U) - | T2J.1}] 

+ 0-j=[\ T2I2) -| T1I1)], 



a 2 = 



At 2 



E 2 + 4t 2 



a 2 +(3 2 = 1. 



(37) 
(38) 



Moreover, only the antisymmetric combination of the 
spin-flip operators 



1 



(39) 



induces a transition to a state with energy E — 0, so that 
the excitation energy is given by Wgpin ~ ~Eo- 

Notice that the exact solution does not display a phase 
transition but remains paramagnetic (and analytic) in 
U/t. On the other hand, in the HF theory, one finds a 
paramagnetic solution below U^.f t jt — 2 and a N eel-type 
ordered solution for U > U^.f t . The latter is clearly non- 
physical and related to the mentioned limitation of mean- 
field in a low-dimensional system. In the GA approxima- 
tion, the electronic correlations are approximated in a 
better way and the range of the paramagnetic solution is 
extended, giving rise to: 

C/ c G iA = 8(V2-l)«3.31 

Since the analytic expressions for the symmetry-broken 
regime become quite lengthy, we restrict the derivation 
below to the paramagnetic case, where the expansion of 
the energy functional is given by 



SE spin 



4C^ 
N 



2^ 8S+SSZ q , 



q—0,7r 

(2 + u)(1-m) 
1 + w 



U/(8t). 



Note that within the HF+RPA approximation, we have 
U s = -U/A. 

The RPA matrices read as: 



.4 



AE + 2U S 

AE + 2U S 



B 



2U S 
2U S 



with the one-particle excitation energies AE = 2£(1 — u 2 ). 

The diagonalization of the eigenvalue problem yields 
two degenerate excitation energies: 



u\ =1 2 = n 2 = AE[AE + 4U S }. 



(40) 



Since the ground state is a singlet, these energies in the 
spin channel coincide with the longitudinal magnetic ex- 
citations computed in the charge channel (see Ref . Il7^ . 
Correspondingly, one has three triplet excitations in to- 
tal, with AS Z = —1,0,1. This indicates that the spin- 
rotation invariance is correctly implemented in our ap- 
proach, a fact that is far from being trivial. It is worth 
noting that in the charge channel an extra assumption 
was needed, namely the antiadiabatic adjustment of the 
double occupancy to the time evolution of the density 
matrix, which was not necessary for the present calcula- 
tion. Therefore, the fact that the spin-rotation symmetry 
is preserved among both independent computations can 
be used as a justification a posteriori of the previous as- 
sumption. 
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FIG. 1: Spin excitation energies for a two-site system com- 
puted with GA+RPA, HF+RPA and the exact result. In the 
inset: J °° Imxt~ {u)duj for the same methods. 



Fig. n shows the magnetic excitation energies for the 
GA+RPA, the HF+RPA and the exact solution. Al- 
though in the exact solution there is no transition at 
finite U and the magnetic excitation energies are always 
finite, in the mean-field approximations there is a tran- 
sition to a symmetry-broken spin-dcnsity-wavc (SDW) 
state, marked by the vanishing of £1 at the respective 
critical values Uf r f/ HF . 

The dynamical transverse susceptibility is given by 



1 ^ (RPA\S+\X)(X\S-\RPA) 



A 



uj — u)\ + le 



(41) 



which is only different from zero for q = tt since 

12 _ — x ( 42 ) 



A F 1 

"£\(RPA\S±\\)\ 2 = ^r6 a 



2VL 



In the inset of Fig.^ we compare the integrated suscepti- 
bility J °° Imxt~{u)duj for the HF+RPA, the GA+RPA 
and the exact result as a function of U ft. The transi- 
tion to the SDW state in the GA and HF approxima- 
tion is characterized by diverging spin-spin correlations 
at q = it, which is due to the vanishing of f2 at Uf r f/ HF . 

The error in both the excitation energies and the inte- 
grated susceptibility is more than halved in the GA+RPA 
with respect to the HF+RPA. This can be traced back 
to the fact that, within the GA, we have a better ground 
state with a delayed instability as a function of U/t. 
Close to the instability one has the worst performance in 
a mean-field plus RPA approach since in this region the 
harmonic approximation to the energy functional breaks 
down. 



B. Paramagnetic regime in infinite dimensions 

As a further application and to get more insight into 
our approximation, we apply the GA+RPA method to an 



infinite-dimensional hypercubic lattice, where the perfor- 
mance is expected to be the best. We consider a partially 
filled system with density n = 1 — S. 

The on-site elements of the density matrix for a param- 
agnetic saddle-point solution are given by plf — j5 aa > , 
so that the matrix Zj of Eq. © reads as 



Z; = 



Z 

z 



(43) 



where, by using the notation introduced by Vollhardt in 
Ref. 0, we have: 



2x 2 



6 2 



1-6 2 



D 



D. 



(44) 
(45) 



For the Gutzwiller approximated energy one obtains 



GA 



E 



zle 



Nz 2 e + NUD, 

1 \ " cr,cr 



(46) 
(47) 



where eo denotes the energy per site of the non- 
interacting system, is the electronic dispersion cor- 
responding to the Gutzwiller Hamiltonian (|19f) and N is 
the number of sites. The minimization of Eq. I|46|) yields 



x i {l-x 2 ) 
x A — S 2 



(1 



U 



(48) 



which, by using Eq. I|45(l . determines the double- 
occupancy parameter D. 

The energy expansion Eq. (|2(J|) in the momentum space 
is given by: 



§E s Pm = ±j2 N « ss i 6S - q 

' *'E (ST+SSZ q + 5S+STI 9 



Nz 

with the following definitions: 

N q = 2e z z" 

k 

6T q ° = $> fe±g + £fe «+" 8 % 



N 4^ 



£ k+qPkk 



(49) 

(50) 
(51) 
(52) 



and the derivatives z' and z" are given in appendix lAl 

Within the RPA approach presented in Sec.[H] one al- 
ways computes all excitation energies, which constitutes 
a suitable procedure for the solutions on finite clusters. 
In infinite system it is usually more convenient to treat 
the RPA problem in terms of a conventional Dyson ap- 
proach. Therefore, we use the well known equivalence 
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between both formulations^ to set up a Dyson equa- 
tion. The interaction kernel which enters the S-matrix 
in the Green's function description can be formally ob- 
tained from Eq. H49[) by substituting the density matrix 
fluctuations by the corresponding operator expressions, 
for instance: 



S q ~ XM'+ 9 ,T Cfe .i-' 
fc 

T q = ^2( £ k + q + e k )cl + ql C ktl . 



Since the energy expansion Eq. (|49|l is a quadratic form 



m 6S± and 8Tf it is useful to define the following matrix 
for the bare time-ordered correlation functions 



jim - i_( (TS+(t)SZ q (0)) (TS+(t)TI q (0)) o 
W ~ N { (7T+(t)S- (0)) (TT+(t)TZ q (0)) 

(53) 

where, the notation (. . . )o indicates that the correlation 
functions are calculated from the excitation spectrum of 
the Gutzwiller Hamiltonian Eq. Ijl9|) and (|29|l and as a 
function of frequency one obtains: 



0/ \_ __Ly^| ^ £ k+£k+q 

W > ~ N V V £ k + £ k+ q (£k + £k+ q ' 2 

k x 



n k+qA {l - n k<i ) n kd (l - n k+qA ) 



LO + £ k+q - £ k +i5 UJ + £ k+q -ek—iS 



r 



(54) 



The RPA series for the spin excitations then corresponds 
to the following Dyson equation: 



with the interaction kernel: 



M„ 



±- 



(55) 



(56) 



As a check of the consistency of our approach, we deter- 
mine the paramagnetic-ferromagnetic and paramagnetic- 
antifcrromagnetic phase boundaries. This can be com- 
pared with previous results within the GA obtained by 
evaluating the vanishing of the corresponding order pa- 
rameter. In case of the ferromagnetic instability we have 
to analyze the limit lim 9 „»o Xg( u — 0) so that the suscep- 
tibility matrix simplifies to 



X° (to = 0) = N(s F ) 



1 2e F 

2e F 4e% 



(57) 



where N(ef) denotes the density of states at the Fermi 
level £p. The inversion of Eq. (|55|l yields as a condition 
for the existence of a pole at u> = and q = 0: 



Det [1 + x °(uj = 0)M ] =1 + ^ = 0, 
with the Landau parameter Fq 



K = N(ef) 



e (2z z'' 



z' 2 ) + 4e F - 



(58) 



(59) 



In the half-filled case (5 = 0) and a symmetric density of 
states (ef = 0) this expression naturally coincides with 
Vollhardt's result [see Eq.(61) in Ref.0]. Fig. [3 displays 
Fq for a Gaussian density of states 



2ttB 



exp 



-—) 

2B 2 J ' 



(60) 
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FIG. 2: Landau parameter Fq as function of U/B for an 
infinite-dimensional hypercubic lattice. The inset shows the 
paramagnetic-ferromagnetic instability line. 



which corresponds to an infinite-dimensional hypercubic 
lattice. In this case the GA becomes exact for the energy 
functional of the Gutzwiller wave function. Due to the 
occurrence of the Brinkman-Rice transition at half filling 
Fq saturates at a value Fq > — 1 for U > 1. Thus, in this 
particular case, there is no second-order paramagnetic- 
ferromagnetic phase transition. The condition Fq = — I 
can be fulfilled in a restricted doping range, i.e., < 5 < 
0.418, and the corresponding instability line is shown in 
the inset of Fig. [2J We find complete agreement of our 
RPA approach with the phase diagram determined by a 
variational approach in Ref. l29l 

In order to investigate the instability toward antiferro- 
magnetism, we study the u> = susceptibility at wave 
vector Q — {%, n, tt, . . . ). The inspection of Eq. I|54|) 
reveals that in the case of a nearest-neighbor hopping 
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FIG. 3: RPA static susceptibility [xq(w = 0)] 11 as function of 
U/B for an infinite dimensional hypercubic lattice. The inset 
shows the paramagnetic-antiferromagnetic instability line. 



netic and antiferromagnetic phases. Since our inten- 
tion is limited to a demonstration of the consistency of 
the GA+RPA approach, we refer the reader to Ref. IH, 
where the antiferromagnetic-ferromagnetic phase bound- 
aries have been determined by comparing the respective 
ground-state energies. 

C. Comparison with exact results 

In the previous subsection we have mainly focused on 
the static limit of our RPA approach. This final part is 
devoted to an analysis of the magnetic properties of the 
GA+RPA method, which is compared to the HF+RPA 
and the exact results on a 4 x 4 Hubbard cluster with 
nearest-neighbor hopping. 



tight-binding band with Ek+ q = —£k only the (1, 1) ma- 
trix element of the bare susceptibility is different from 
zero: 



1 



87TB 



-Et 



14 

2 B 2 



(61) 



where E\ (x) denotes the exponential integral. 30 The RPA 
series of Eq. 1|55JI then leads to 



[XQ(w = 0)] u = 
Nq = 



[x° (" = o)] n 



l + iV Q [x° Q (u = 0) 
e {2z z" - (z 1 ) 2 }. 

0)] n for^ 



(62) 



li 



We show the behavior of [xq(^ 
Fig.El 

Due to the complete nesting, the bare susceptibility 
[xq( w = 0)] diverges for 5 = 0. Hence in this case the 
singularities of [xq{& = 0)]ii are determined by the zeros 
of the interaction kernel Nq which naturally vanishes for 
U /B = but also at the Brinkman-Rice transition where 
zq — > 0. The latter, however, is irrelevant since it occurs 
in the antiferromagnetic phase. The pole at U/B = in- 
dicates that the instability toward antiferromagnetism at 
half filling occurs at arbitrarily small interaction also in 
infinite dimensions. For finite 5 the bare magnetic suscep- 
tibility is finite and consequently the pole of \Q (k> = 0) is 
due to the vanishing of the RPA denominator in Eq. (|62() . 
It turns out that the static magnetic susceptibility has ex- 
actly one pole in the range < 6 < 0.117, two poles in the 
range 0.117 < S < 0.2048 and no pole for S > 0.2048. For 
completeness, Fig.|3]also displays xq( w = 0) for 8 = 0.25, 
where there is a small enhancement for those values of 
U/B where the instability occurred for smaller 5. The 
inset of Fig.[3]shows the antiferromagnetic-paramagnetic 
instability line constructed from the poles of xq( w = 0). 
Again we find complete agreement with the variational 
approach of Ref. |22. Note that one should also deter- 
mine the first-order boundaries between the ferromag- 



1. Half -filled system 

We start with the half-filled system with 8 spin-up and 
8 spin-down particles. The ground state Slater determi- 
nant for the GA and the HF approximation corresponds 
to a SDW, which breaks the spin-rotational symmetry of 
the Hamiltonian. As a consequence the transverse mag- 
netic excitations contain zero-energy Goldstone modes at 
wave vector Q = (ir,ir). To avoid numerical instabilities, 
we have added a small perturbation to the Hamiltonian 



V = 



(63) 



with a ~ 10~ 4 £, which shifts the Goldstone modes to 
small but finite energies (~ a). In the exact solution an 
analogue pole appears at small but non-zero frequency 
(u>/t w 0.145) due to the finiteness of the cluster. In 
the thermodynamic limit long-range order is recovered 3 ^ 
and a Goldstone mode will appear as in the mean-field 
solution with a weight related to the order parameter. 
Here, we are interested in the finite-frequency behavior 
and, therefore, we exclude the exact and approximate 
"Goldstone-like" poles from the comparison and restrict 
ourself to the finite- frequency (triplet) excitations, which, 
for the chosen value of a, do not sensitively depend on 
the anisotropy field Eq. (j^SJ- 

Fig.0]shows the magnetic excitation energies as a func- 
tion of U/t evaluated within the GA+RPA, the HF+RPA 
and the exact diagonalization. Note that the 4x4 system 
has a further accidental symmetry, which causes degen- 
eracy between the q = (tt/2, 7r/2) and q = (tt, 0) excita- 
tions. Furthermore, the SDW ground state of the GA and 
HF solution leads to the doubling of the Brillouin zone 
(see inset of Fig.0J) so that, besides the antiferromagnetic 
wave vector Q, only q = (tt/2, 0) and q = (tt, 0) corre- 
spond to independent excitations. On the other hand, 
on the 4x4 lattice, we have that the exact energies at 
q = (7r/2,0) and q = (n/2,Tr) are slightly different. 

The small-U behavior of the lowest excitation energy 
in Fig. 0| can be well understood from the SDW picture. 
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FIG. 4: Magnetic excitations at q — (it/2, 0) and q — (tt, 0) as 
a function of U/t for a half-filled 4x4 cluster: GA+RPA (solid 
line), HF+RPA (dashed line), and exact diagonalization (full 
circles and full squares). The exact diagonalization results 
for the excitation at q — (n/2,n) are also reported (empty 
squares). In the inset: the q-point mesh of the 4x4 cluster 
and the dashed square indicates the doubled Brillouin zone. 
Shaded points indicates the important wave vectors of the 
magnetic excitations. 



Within this approximation, the band structure in the re- 
duced Brillouin zone is given by E q = ± . / e 1 + A 2 , with 



e q = — 2t[co$(q x ) + cos(q y )] and A denotes the SDW gap. 
Since we study a half-filled system, all states with E q < 
are occupied. Consider first the q = (ir, 0) excitation 
which can be attributed to a spin-flip transition from 
qx = (-71-/2, ±7r/2) to q 2 = (7r/2,±7r/2) so that the ex- 
citation energy is given by lu = E qi — E q2 — 2 A. The 
SDW gap in the HF approximation is related to the on- 
site magnetization A HF = 2f7|S' z |, whereas within the 
KR formulation of the GA it is determined by the dif- 
ference in the local spin-dependent Lagrange multipliers 
A GA = A| — A j,. Since in the limit U — > the GA reduces 
to the HF approximation, both excitation energies coin- 
cide in this regime and also agree with the exact result. 
On the other hand, for U/t > 1 , where RPA corrections 
become important, it can be seen from Fig. Q]that the 
GA+RPA is in much better agreement with exact diag- 
onalization than the corresponding HF+RPA result. As 
a consequence, the GA+RPA gives a quite accurate de- 
scription of the crossover (at U/t » 6) from the SDW 
regime, where a gap proportional to U opens along the 
Fermi surface, to the Heisenberg regime, where there are 
low-energy magnetic excitations with energy scale t 2 /U . 

For the higher energy triplet excitation at q = (tt/2, 0), 
the GA+RPA yields energies which are slightly lower 
than the exact result. However, whereas the discrep- 
ancy for the GA+RPA at U/t = 6 is around 10%, the 
HF+RPA deviates by almost 20% from the exact diago- 
nalization result. 



FIG. 5: Local magnetic susceptibility x{ u ) f° r the half- 
filled 4x4 cluster calculated within exact diagonalization, 
GA+RPA and HF+RPA for U/t = 4. The HF+RPA curve 
has been shifted for convenience. The two arrows indicate 
the energy of the lowest Q — (tt,tt) excitation at u/t ~ 0.145 
(exact diagonalization) and u/t = (RPA Goldstone mode). 



In order to have information on the accuracy of the 
GA+RPA for finite frequencies, we report in Fig. [5] the 
local magnetic susceptibility 

XM = E E l(*m|S+|*o>| 2 <^ - (E m - E )), (64) 

q m>0 

for the GA+RPA and the HF+RPA approximations and 
the exact diagonalization for U/t = 4. The 5-functions in 
Eq. (|64|l have been replaced by Lorentzians with width 
O.lt. 

The two lowest-energy excitations are quite accurate 
within the GA+RPA approach except for a moderate 
overestimation of the intensity which becomes much 
worse in HF+RPA. Interestingly, also the high-energy 
spin fluctuations are in the correct frequency range. 

It should be noted that the spectral weight is con- 
strained by the following sum rule: 



duux(v) = --{T)ga, 
o 1 



(65) 



where {T)qa is the average value of the kinetic energy 
in the GA. The sum rule Eq. l|f)K|) relates the RPA cor- 
relation function x( w ) to the kinetic energy computed 
within the GA. A similar sum rule is valid in HF+RPA 
with the kinetic energy computed in HF. In Ref. Il6l we 
have already demonstrated that the GA kinetic energy 
is in remarkable agreement with the exact result over a 
large doping range which in the present context gives ad- 
ditional support to the GA+RPA approach also in the 
magnetic sector. On the other hand, the HF approxima- 
tion is of inferior quality in describing excitation energies 
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FIG. 6: Cumulative sum of the first moment of x(^) f° r the 
exact result, GA+RPA, and HF+RPA. Data are for a half- 
filled 4x4 cluster and U/t = 4. Inset: a detail of the low- 
energy part. 

and the total kinetic energy. Therefore it is not surpris- 
ing that also spectral weights perform much worse than 
in the GA+RPA approach. 

Finally, Fig. [fj] displays the frequency evolution of the 
first moment of x(u;): 

AfV) = / (1CjCjx{Cj), (66) 
Jo 

for the same parameters as in Fig. EI Note that M 1 (w) 
contains the contribution from the lowest Q = (tt,tt) 
excitations (i.e., the Goldstone modes in the mean- field 
plus RPA), which appear as the offsets at small energies. 
From the inset of Fig. it can be seen that, especially 
at low frequencies, the GA+RPA approach provides a 
much better approximation for the exact M 1 (w) than 
HF+RPA, both with respect to the poles and intensities. 
Both the HF+RPA and the GA+RPA approximate the 
incoherent part of the exact spectrum (i.e., for u>/t > 2) 
by a rather small set of excitations. However, the corre- 
sponding step-like evolution of the first moment of x( w ) 
is quite close to the exact result. 



2. Doped system 

We finally investigate a 4 x 4 cluster with 5 spin-up 
and 5 spin-down particles, corresponding to a closed-shell 
configuration. The HF solution undergoes a magnetic in- 
stability with wave vector q = (tv/2, tt) at U cr it/t w 4.365, 
marked by the softening of the corresponding excitation 
at this critical value. For U > U cr it the HF+RPA spec- 
trum has a Goldstone mode and in general the perfor- 
mance is very poor consistently with the fact that a bro- 
ken symmetry state is not expected even in the thermo- 



4 



\\ 
'■ v\ 
\ \\ 

\ \\ q=(u,ic) 
i \ ^\ 
\ ^\ 

N. ^^^^ 


Y\ 

- V \ — exact 

\ \ -- GA+RPA 
. \ ''■ \ HF+RPA . 

\ \ 
\ \ 

\ N 
\ X 


\v 


\ \. 
\ 

\ — 


\ X 

\ \^ q=(Jt,0) 

" \ Yv 

\ \\ 
■•• ^\ 

^ x^ 
'■ ^ 




\ — — . 

\ s x q=(3i/2,3t) 

'. N. - 


\ • — 

C ' Mi 



2 4 6 8 1012141618 2 4 6 8 10 12 14161820 
U/t U/t 

FIG. 7: Magnetic excitations for wave vectors q = (n, n) 
and q = (n/2,n) (left panel) and for q — (7r,0), q = (tt/2,0) 
(right panel) as a function of U /t for the exact diagonalization 
(solid line), GA+RPA (dashed line), and HF+RPA (dotted 
line). The HF+RPA are shown for U < U cr it- 



dynamic limit, therefore, we restrict to the more physi- 
cal paramagnetic solution. In the GA case the param- 
agnetic solution can be stabilized for much larger val- 
ues of U/t providing a reasonable starting point in a 
broader parameter range just as in the two-site case. The 
GA+RPA approach captures the behavior of the exact 
solution (namely the softening of triplet excitations) at 
least in a qualitative way, although quantitative devia- 
tions increase with increasing U/t (see Fig.[7J|. 

In Fig. we compare the local susceptibility of the 
HF+RPA, the GA+RPA, and the exact diagonalization 
for U/t — 4, i.e., for values of U/t smaller than the mag- 
netic instability. The GA+RPA not only gives a rather 
good estimate to the lowest excitation energy but in addi- 
tion provides a good approximation for the corresponding 
intensity. Note that, since for the given value of U/t the 
HF solution is already close to the q = (ir/2,ir) insta- 
bility, we observe a strong softening of the lowest energy 
excitation resulting in a significantly enhanced oscillator 
strength. 

Finally, we have also evaluated the cumulative integral 
of the first moment of x( w ) f° r the GA+RPA and exact 
diagonalization, which is shown in the inset of Fig. [S] 
Due to the sum rule Eq. (|65|l and the excellent kinetic 
energy approximation of the GA, the integrated weight 
of the GA+RPA and the exact diagonalization are in 
excellent agreement. Moreover, we again observe that the 
GA+RPA provides a rather good step- like approximation 
to the exact evolution of the spectral weight as a function 
of the frequency. 
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FIG. 8: Local magnetic susceptibility x(lo) for the 4x4 clus- 
ter with 10 particles calculated within exact diagonalization, 
GA+RPA and HF+RPA for U/t = 4. The HF+RPA curve 
has been shifted for convenience. The inset shows the cumu- 
lative sum of the first moment of xi. 1 - ) f° r the exact result 
and the GA+RPA. 



IV. CONCLUSION 

In this paper, we have presented a detailed investi- 
gation of the quality of the GA+RPA approach for the 
computation of the magnetic excitations. The present 
computation is complementary to the previous compu- 
tation in the charge sectoriiSiii An unexpected outcome 
of the present work is a further justification of the antia- 
diabatic assumption for the time-evolution of the double 
occupancy, which was needed in the charge but not in 
the magnetic channel. The fact that the two calculations 
for the spin and the charge sector give the correct de- 
generacy of the excitation spectrum for a spin-rotational 
invariant system (the two-site Hubbard model) clearly 
indicates that such an assumption was indeed correct, 
i.e., other possibilities like to keep the double occupancy 
fixed at the stationary value (rather than to follow the 



time evolution of the density matrix) would had lead to 
an unphysical breaking of spin-rotational symmetry. 

The present formalism is based on a Gutzwiller-type 
energy functional, which can be either obtained from 
the spin-rotational invariant KR slave-boson scheme 2 ^ 
or alternatively from the standard GA with spin-rotated 
Slater determinants. 24 In our approach, due to the fact 
that all bosonic fields have been already eliminated from 
the saddle-point energy functional, it turns out that the 
evaluation of RPA fluctuations around the GA solution is 
significantly simplified. In the present paper, we have re- 
stricted the calculations of the magnetic excitations to 
small Hubbard clusters in order to compare with ex- 
act diagonalization results. The better performance of 
GA+RPA with respect to HF+RPA has been demon- 
strated for both excitation energies and the correspond- 
ing intensities. However, compared to numerical meth- 
ods^ our approach can be pushed to much larger sys- 
tems. In particular, it is suitable for the evaluation 
of magnetic excitations around inhomogeneous solutions 
of Hubbard-type models, where it is constrained to the 
same size limitations than the unrestricted HF+RPA ap- 
proach. This is interesting in connection with the mag- 
netic susceptibility in nickelates and high-T c cuprates, 
which are both characterized by the presence of strong 
electronic correlations and inhomogeneous charge distri- 
butions in some part of the phase diagram. Work in this 
direction is in progress. 
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APPENDIX A: DERIVATIVES OF THE 
HOPPING FACTOR 

The derivatives appearing in Eqs. H21|) and lll'L'l) are 

jiven by: 
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In case of a homogeneous, paramagnetic saddle point these expressions simplify to 
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where zq denote the elements of the z-matrix defined in where 
Eq. igU). 
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APPENDIX B: DERIVATION OF THE RPA 

MATRICES dZi U dz 



' 9(Pjj,n) d {pjj,ir) 



In order to give explicit expressions for the RPA ma- 
ices A and £? as defined in Eq. (|30|) one has firs 
diagonalize the Gutzwiller Hamiltonian Eq. 119(1 via 

^2*i(v,o)a Vtt „ (Bl) 



trices A and £> as defined in Eq. (|3"0|) one has first to Tp£ p , h , = 5 aa >5 a ]- V^tjj ^ — T)^j(^ I) 

• • - ' y "Pjj,U 

x 0? n <fi(p' T)$i(V 4) + < u $ i (p' 4)) 

where i/ refers to either particle (p) or hole (h) states. / / Q / f 

Inserting this transformation in the expansion Eq. jgQ) x { z i,n®j(P l) $ *0' T) + \u^iip l)®j(h T)) , 

leads to the following expressions for A and B: . . \ it . , . , _ ,, , , 

Klp'h' = ^,-^^TS*<i^r ±L *j(pt)*i(fti) 



8„„>^N ij $ i fa)$ i (h,-a)$ j (j3 , cT)$j(h',-o-) 



+ K ph,p'h> + U p >h>.ph' 3 

« ( B3 ) 
x ^ (p', -ff)^ (/iV)) + T;-;,,, + Tp,^, ph , (B2) 
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